function makeRunData_consumerAdoptions(runId)

addpath('../general_funcs');

%%% Read config for run
config = getConfig_consumerAdoptions(runId);

if config.modelId == 3
	makeRunData_consumerAdoptions_ctsTime(runId);
	return;
end

dataId = config.dataId;

spatialMatricesFolder = config.spatialMatricesFolder;
runDataFile           = config.runDataFile;
dynPredInputFile      = config.dynPredInputFile;

%%%%% Get covariates (including FEs)
covariates = config.covariates;

spatialMatrices = covariates.spatialMatrices; % (k1-k2) covariates

myX_tk1_names   = covariates.myX_tk1_names;
myX_tk1_transfs = covariates.myX_tk1_transfs;
		
myX_tk2_names   = covariates.myX_tk2_names;
myX_tk2_transfs = covariates.myX_tk2_transfs;
		
myZ_k1_names   = covariates.myZ_k1_names;
myZ_k1_transfs = covariates.myZ_k1_transfs;
		
myZ_k2_names   = covariates.myZ_k2_names;
myZ_k2_transfs = covariates.myZ_k2_transfs;
		
myAnames   = covariates.myAnames;
myAtransfs = covariates.myAtransfs;
		
myZ_k1_geo_FE_Names = covariates.myZ_k1_geo_FE_Names;
myZ_k2_geo_FE_Names = covariates.myZ_k2_geo_FE_Names;
time_FE_Names = covariates.time_FE_Names;


%%% Load files
data1 = loadDataFiles(dataId);

NewConsumers        = data1.NewConsumers;        % (T*K1*K2) x 1 (sparse)
NewProviders         = data1.NewProviders;         % T x K
DeletedProviders     = data1.DeletedProviders;     % T x K
laggedConsumersOrig = data1.laggedConsumersOrig; % T x K1
laggedConsumersDest = data1.laggedConsumersDest; % T x K2
laggedProviders      = data1.laggedProviders;      % T x K
population        = data1.population;        % K x 1
surfaces          = data1.surfaces;          % K x 1
incomingCommuters = data1.incomingCommuters; % K x 1
touristBeds       = data1.touristBeds;       % K x 1
periodLabels = data1.periodLabels;
geoLabels0   = data1.geoLabels;
A     = data1.A;
Z     = data1.Z;
X     = data1.X;
Anames     = data1.Anames;
Znames     = data1.Znames;
Xnames     = data1.Xnames;
surfaceName         = data1.surfaceName;
geoAggLevelNames    = data1.geoAggLevelNames;
geoAggLevelLabels   = data1.geoAggLevelLabels;
timeAggLevelNames   = data1.timeAggLevelNames;
timeAggLevelLabels = data1.timeAggLevelLabels;

% Read dimensions
[T,K] = size(NewProviders);
K1 = K;
K2 = K;

%%% Make Y (=NewConsumers)
Y = NewConsumers;
clear NewConsumers;

%%% Make popAtRisk and logLambdaOffset
marketSize = Z(:,find(strcmp(Znames, 'NumPersons19orOlder'))); % K x 1
initialSeed = laggedConsumersOrig(1,:)'; % K1 x 1 
popAtRisk = marketSize' - laggedConsumersOrig; % T x K1
logLambdaOffset_dim1 = log(popAtRisk); % T x K1

% Reshape installedBases
laggedProviders      = reshape(laggedProviders, [T*K2 1]);      % (T*K2) x 1
laggedConsumersOrig = reshape(laggedConsumersOrig, [T*K1 1]); % (T*K1) x 1
laggedConsumersDest = reshape(laggedConsumersDest, [T*K2 1]); % (T*K2) x 1

%%%%%%%%%% TRICK TO HANDLE CASE WHERE laggedProviders=0 %%%%%%%%%%
% Trick: add a "penalty" to logLambdaOffset_{tk2} where tk2 is "bad"
%	--> with this penalty, logLambda will be very negative, which means that lambda will be effectively
% (computationally speaking exactly) equal to zero. Yet, V is will not be -Inf and that is a good thing, because
% otherwise, when doing Y.*V we have a problem (Y equal to 0, V equal to -Inf, the product gives NaN).
% Similar issues arise with the computation of the gradient and of the hessian, which involve the computation of X.*(Y-lambda)
% and X.*X.*lambda. In those cases, X cannot be equal to -Inf either (which would happen if we take beta*log(laggedConsumersOrig)),
% so the trick is to set X to some other value. Any non-infinite value for X will do, given that (Y-lambda) = (0-0) = 0 for the grad,
% and lambda=0 for the hessian.
is_bad_tk2  = laggedProviders <= 0; % (T*K2) x 1
logLambdaOffset_dim2   = -2e10 * reshape(is_bad_tk2, [T 1 K2]); % T x 1 x K2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dimType = 3;

NumTK2 = size(laggedProviders,1);

% Get myX_tk1 and myX_tk1_names
[Xflag, X_idces] = ismember(myX_tk1_names, Xnames);
if any(~Xflag); error('Could not find all Xnames.'); end;
if any(~strcmp(Xnames(X_idces), myX_tk1_names)); error('Mistmatch'); end;
myX_tk1         = X(:,:, X_idces);				% T x K x NumX
myX_tk1         = reshape(myX_tk1, T*K, size(myX_tk1,3));	% (T*K) x NumX
[myX_tk1, myX_tk1_names] = transformMatrix(myX_tk1, myX_tk1_transfs, myX_tk1_names, reshape(repmat(surfaces', T, 1), T*K, 1),surfaceName);
	% (T*K) x NumX and cell(1,NumX) of strings

% Get myX_tk2 and myX_tk2_names
[Xflag, X_idces] = ismember(myX_tk2_names, Xnames);
if any(~Xflag); error('Could not find all Xnames.'); end;
if any(~strcmp(Xnames(X_idces), myX_tk2_names)); error('Mistmatch'); end;
myX_tk2         = X(:,:, X_idces);				% T x K x NumX
myX_tk2         = reshape(myX_tk2, T*K, size(myX_tk2,3));	% (T*K) x NumX
[myX_tk2, myX_tk2_names] = transformMatrix(myX_tk2, myX_tk2_transfs, myX_tk2_names, reshape(repmat(surfaces', T, 1), T*K, 1),surfaceName);
	% (T*K) x NumX and cell(1,NumX) of strings

% Get myZ_k1 and myZ_k1_names
[Zflag, Z_idces] = ismember(myZ_k1_names, Znames);
if any(~Zflag); error('Could not find all Znames.'); end;
if any(~strcmp(Znames(Z_idces), myZ_k1_names)); error('Mistmatch'); end;
myZ_k1         = Z(:, Z_idces);				% K x NumZ
[myZ_k1, myZ_k1_names] = transformMatrix(myZ_k1, myZ_k1_transfs, myZ_k1_names, surfaces, surfaceName);
	% K x NumZ and cell(1,NumZ) of strings

% Get myZ_k2 and myZ_k2_names
[Zflag, Z_idces] = ismember(myZ_k2_names, Znames);
if any(~Zflag); error('Could not find all Znames.'); end;
if any(~strcmp(Znames(Z_idces), myZ_k2_names)); error('Mistmatch'); end;
myZ_k2         = Z(:, Z_idces);				% K x NumZ
[myZ_k2, myZ_k2_names] = transformMatrix(myZ_k2, myZ_k2_transfs, myZ_k2_names, surfaces, surfaceName);
	% K x NumZ and cell(1,NumZ) of strings

% Get myA and myAnames
[Aflag, A_idces] = ismember(myAnames, Anames);
if any(~Aflag); error('Could not find all Anames.'); end;
if any(~strcmp(Anames(A_idces), myAnames)); error('Mistmatch'); end;
myA         = A(:, A_idces);				% T x NumA
	% T x NumA and cell(1,NumA) of strings

% Add intercept to myA
myA = [ones(T,1) myA];
myAnames = [{'Intercept'} myAnames];
myAtransfs = [-1 myAtransfs];

%% Make time_FEs
[flag, idxes] = ismember(time_FE_Names, Anames);
if any(~flag); error('Could not find all time FE names.'); end;
time_FEs = A(:,idxes);
Num_time_FE_sets = length(time_FE_Names);
Num_time_FE_vals = zeros(1,Num_time_FE_sets);
time_FE_labels = cell(1,Num_time_FE_sets);
for aa = 1:Num_time_FE_sets
	time_FEs_aa = time_FEs(:,aa);
	time_FE_Name_aa = time_FE_Names{aa};
	
	unq_time_FEs_aa = sort(unique(time_FEs_aa));
	time_FE_labels_aa = strcat(time_FE_Name_aa, sprintfc('_is_%d', unq_time_FEs_aa));
	[flag,newTime_FEs_aa] = ismember(time_FEs_aa, unq_time_FEs_aa);
	
	time_FEs(:,aa) = newTime_FEs_aa;
	Num_time_FE_vals(aa) = length(unq_time_FEs_aa);
	time_FE_labels{aa} = time_FE_labels_aa;
end

%% Make  myZ_k1_geo_FEs + get myZ_k1_geo_FE_labels +  myZ_k1_Num_geo_FE_vals
[flag, idxes] = ismember(myZ_k1_geo_FE_Names, Znames);
if any(~flag); error('Could not find all geo FE names.'); end;
myZ_k1_geo_FEs = Z(:,idxes);
Num_geo_FE_sets = length(myZ_k1_geo_FE_Names);
myZ_k1_Num_geo_FE_vals = zeros(1,Num_geo_FE_sets);
myZ_k1_geo_FE_labels = cell(1,Num_geo_FE_sets);
for aa = 1:Num_geo_FE_sets
	geo_FE_Name_aa = myZ_k1_geo_FE_Names{aa};
	[flag,geoidx] = ismember(geo_FE_Name_aa, geoAggLevelNames);
	if ~flag; error('Could not find geo_FE_name among geoAggLevelNames'); end;
	geo_FE_labels_aa = geoAggLevelLabels{geoidx};
	myZ_k1_Num_geo_FE_vals(aa) = length(geo_FE_labels_aa);
	myZ_k1_geo_FE_labels{aa} = geo_FE_labels_aa;
end

%% Make  myZ_k2_geo_FEs + get myZ_k2_geo_FE_labels +  myZ_k2_Num_geo_FE_vals
[flag, idxes] = ismember(myZ_k2_geo_FE_Names, Znames);
if any(~flag); error('Could not find all geo FE names.'); end;
myZ_k2_geo_FEs = Z(:,idxes);
Num_geo_FE_sets = length(myZ_k2_geo_FE_Names);
myZ_k2_Num_geo_FE_vals = zeros(1,Num_geo_FE_sets);
myZ_k2_geo_FE_labels = cell(1,Num_geo_FE_sets);
for aa = 1:Num_geo_FE_sets
	geo_FE_Name_aa = myZ_k2_geo_FE_Names{aa};
	[flag,geoidx] = ismember(geo_FE_Name_aa, geoAggLevelNames);
	if ~flag; error('Could not find geo_FE_name among geoAggLevelNames'); end;
	geo_FE_labels_aa = geoAggLevelLabels{geoidx};
	myZ_k2_Num_geo_FE_vals(aa) = length(geo_FE_labels_aa);
	myZ_k2_geo_FE_labels{aa} = geo_FE_labels_aa;
end


%%% Add treatment for initial installed base %%%
if config.initialBaseTreatment
	[myZ_k1, myZ_k1_names] = addInitialBaseCovariate(initialSeed, myZ_k1, myZ_k1_names, 'Consumers');
end


%%%%% Load spatial matrices
X_k1k2       = zeros(K1*K2,0); % (K1*K2) x NumX_k1k2
X_k1k2_names = {}; % cell(NumX_k1k2)
NumSpatialMatrices = length(spatialMatrices);
for ss = 1:NumSpatialMatrices
	matName = spatialMatrices{ss};
	if strcmp(matName, 'sameCanton')
		W = eye(K); % K x K
	else
		load(sprintf('%s/%s.mat', spatialMatricesFolder, matName), 'W', 'geoLabels');
		if exist('geoLabels')
			if ~isequal(geoLabels, geoLabels0); error('Inconsistency in geoLabels!'); end;
			clear geoLabels;
		end
		W = full(W); % may be needed for sparse matrices (e.g. contiguity matrix)
		if strcmp(matName, 'Weight_TouristTrip') % expressed in million trips
			W = 1e6*W;
		end
		
		if strcmp(matName, 'MobilitePro') % set diagonal equal to zero (otherwise it's not a flow and might even be confused with unemployment)
			W(get_diagonal_idxes(size(W,1))) = 0;
		end
		
		if ~ismember(matName, {'Weight_Contiguity', 'sameDEP'})
			% Normalize W by population
			population = reshape(population, [K1 1]);
			W = W./population;
		end
	end
	W = reshape(W, [K1*K2 1]); % (K1*K2) x 1
	X_k1k2 = [X_k1k2 W];
	X_k1k2_names = [X_k1k2_names {matName}];
	clear W matName;	
end
% need to be careful about scaling (avoid numerical problems)
[X_k1k2, X_k1k2_names] = transformMatrix(X_k1k2, zeros(1,length(spatialMatrices)), X_k1k2_names);




%%%% Make Xdynamic_staticInputs and Xparts_dynamic (Xparts that depend on installed base)
Xparts_static{2}.X            = myX_tk1;       % (T*K1) x NumXstatic
Xparts_static{2}.Xnames       = myX_tk1_names; % cell(1,NumXstatic) of strings
Xparts_static{2}.X_FEs        = zeros(T*K1,0); % (T*K1) x NumX_FEs
Xparts_static{2}.X_FE_labels  = {};            % cell(1, NumX_FEs)
Xparts_static{2}.NumX_FE_vals = zeros(1,0);  % 1 x NumX_FEs

Xparts_static{3}.X            = myX_tk2;         % (T*K2) x NumXstatic
Xparts_static{3}.Xnames       = myX_tk2_names;   % cell(1,NumXstatic) of strings
Xparts_static{3}.X_FEs        = zeros(NumTK2,0); % (NumTK2) x NumX_FEs
Xparts_static{3}.X_FE_labels  = {};              % cell(1, NumX_FEs)
Xparts_static{3}.NumX_FE_vals = zeros(1,0);    % 1 x NumX_FEs

Xdynamic_staticInputs.Xparts_static = Xparts_static;
if config.laggedBaseTransf.code == 1
	Xdynamic_staticInputs.extraArgs     = {};
end
if config.laggedBaseTransf.code == 2
	if strcmp(config.laggedBaseTransf.denom, 'Area')
		surface_rep = reshape(repmat(surfaces', [T, 1]), [T*K 1]); % (T*K) x 1
		Xdynamic_staticInputs.extraArgs = {log(surface_rep)};
	end
end

installedBases.laggedProviders      = laggedProviders;
installedBases.laggedConsumersOrig = laggedConsumersOrig;
installedBases.laggedConsumersDest = laggedConsumersDest;
Xparts_dynamic = make_Xdynamic_consumerAdoptions(config, Xdynamic_staticInputs, installedBases, true);


%%%%% Make Xparts and dims %%%%%
Xparts = {};
dims1 = {};
Xpart_2_NumX = [];
Xpart_2_NumX_FEs = [];
Xpart_2_Num_FE_vals = {};
NumFEvals2Keep = {};

pidx = 0;

%%% K1-K2 part: spatial matrices
X            = X_k1k2;            % (K1*K2) x NumX
Xnames       = X_k1k2_names;      % 1 x NumX
X_FEs        = zeros(K1*K2, 0);   % (K1*K2) x NumX_FEs
X_FE_labels  = {};                % cell(1, NumX_FEs)
NumX_FE_vals = zeros(1,0);        % 1 x NumX_FEs

pidx = pidx+1;
obj.X                     = X;
obj.Xnames                = Xnames;
obj.X_FEs                 = X_FEs;
obj.X_FE_labels           = X_FE_labels;
obj.NumX_FE_vals          = NumX_FE_vals;
Xparts{pidx}              = obj;
dims1{pidx}               = K1*K2;
Xpart_2_NumX(pidx)        = size(X,2);
Xpart_2_NumX_FEs(pidx)    = size(X_FEs,2);
Xpart_2_Num_FE_vals{pidx} = NumX_FE_vals;
NumFEvals2Keep{pidx}      = sum(NumX_FE_vals) - length(NumX_FE_vals);

clear X X_FEs NumX_FE_vals obj;


%%% T-K1 part:
pidx = pidx+1;
obj = Xparts_dynamic{2};
Xparts{pidx}              = obj;
dims1{pidx}               = T*K1;
Xpart_2_NumX(pidx)        = size(obj.X,2);
Xpart_2_NumX_FEs(pidx)    = size(obj.X_FEs,2);
Xpart_2_Num_FE_vals{pidx} = obj.NumX_FE_vals;
NumFEvals2Keep{pidx}      = sum(obj.NumX_FE_vals) - length(obj.NumX_FE_vals);

clear obj;


%%% T-K2 part:
pidx = pidx+1;
obj = Xparts_dynamic{3};
Xparts{pidx}              = obj;
dims1{pidx}               = NumTK2;
Xpart_2_NumX(pidx)        = size(obj.X,2);
Xpart_2_NumX_FEs(pidx)    = size(obj.X_FEs,2);
Xpart_2_Num_FE_vals{pidx} = obj.NumX_FE_vals;
NumFEvals2Keep{pidx}      = sum(obj.NumX_FE_vals) - length(obj.NumX_FE_vals);

clear obj;


%%% T part:
X      = myA;      % T x NumX
Xnames = myAnames; % 1 x NumX
X_FEs        = time_FEs;         % T x NumX_FEs
X_FE_labels  = time_FE_labels;   % cell(1, NumX_FEs)
NumX_FE_vals = Num_time_FE_vals; % 1 x NumX_FEs

pidx = pidx+1;
obj.X                     = X;
obj.Xnames                = Xnames;
obj.X_FEs                 = X_FEs;
obj.X_FE_labels           = X_FE_labels;
obj.NumX_FE_vals          = NumX_FE_vals;
Xparts{pidx}              = obj;
dims1{pidx}               = T;
Xpart_2_NumX(pidx)        = size(X,2);
Xpart_2_NumX_FEs(pidx)    = size(X_FEs,2);
Xpart_2_Num_FE_vals{pidx} = NumX_FE_vals;
NumFEvals2Keep{pidx}      = sum(NumX_FE_vals) - length(NumX_FE_vals);

clear X X_FEs NumX_FE_vals obj;


%%% K1 part:
X      = myZ_k1;        % K1 x NumX
Xnames = myZ_k1_names;  % 1 x NumX
X_FEs        = myZ_k1_geo_FEs;  % K1 x NumX_FEs
X_FE_labels  = myZ_k1_geo_FE_labels;   % cell(1, NumX_FEs)
NumX_FE_vals = myZ_k1_Num_geo_FE_vals; % 1 x NumX_FEs

pidx = pidx+1;
obj.X                     = X;
obj.Xnames                = Xnames;
obj.X_FEs                 = X_FEs;
obj.X_FE_labels           = X_FE_labels;
obj.NumX_FE_vals          = NumX_FE_vals;
Xparts{pidx}              = obj;
dims1{pidx}               = K1;
Xpart_2_NumX(pidx)        = size(X,2);
Xpart_2_NumX_FEs(pidx)    = size(X_FEs,2);
Xpart_2_Num_FE_vals{pidx} = NumX_FE_vals;
NumFEvals2Keep{pidx}      = sum(NumX_FE_vals) - length(NumX_FE_vals);

clear X X_FEs NumX_FE_vals obj;



%%% K2 part:
X      = myZ_k2;   % K2 x NumX
Xnames = myZ_k2_names;  % 1 x NumX
X_FEs        = myZ_k2_geo_FEs;  % K2 x NumX_FEs
X_FE_labels  = myZ_k2_geo_FE_labels;   % cell(1, NumX_FEs)
NumX_FE_vals = myZ_k2_Num_geo_FE_vals; % 1 x NumX_FEs

pidx = pidx+1;
obj.X                     = X;
obj.Xnames                = Xnames;
obj.X_FEs                 = X_FEs;
obj.X_FE_labels           = X_FE_labels;
obj.NumX_FE_vals          = NumX_FE_vals;
Xparts{pidx}              = obj;
dims1{pidx}               = K2;
Xpart_2_NumX(pidx)        = size(X,2);
Xpart_2_NumX_FEs(pidx)    = size(X_FEs,2);
Xpart_2_Num_FE_vals{pidx} = NumX_FE_vals;
NumFEvals2Keep{pidx}      = sum(NumX_FE_vals) - length(NumX_FE_vals);

clear X X_FEs NumX_FE_vals obj;


%%%%% Make data object
% Calculate NumParams
NumParts = pidx;
NumParams = 0;
for ii = 1:NumParts
	NumParams = NumParams + size(Xparts{ii}.X,2);
	for ff = 1:size(Xparts{ii}.X_FEs,2)
		NumX_FE_vals_if = Xparts{ii}.NumX_FE_vals(ff);
		NumParams = NumParams + NumX_FE_vals_if - 1;
	end
end

% Make dims
dims.dimType  = dimType;
dims.NumParts = NumParts;
dims.NumObs = size(Y,1);
dims.NumParams      = NumParams;
dims.dims1               = dims1;
dims.Xpart_2_NumX        = Xpart_2_NumX;
dims.Xpart_2_NumX_FEs    = Xpart_2_NumX_FEs;
dims.Xpart_2_Num_FE_vals = Xpart_2_Num_FE_vals;
dims.NumFEvals2Keep      = NumFEvals2Keep;

% Make LL_offset
myidxes = find(Y>0);
LL_offset = -sum(gammaln(1+Y(myidxes)));
clear myidxes;

% Make data
clear data; % (Clear previous data object)
data.Y               = Y;
data.logLambdaOffset_dim1 = logLambdaOffset_dim1;
data.logLambdaOffset_dim2 = logLambdaOffset_dim2;
data.LL_offset       = LL_offset;
data.Xparts          = Xparts;
data.dims            = dims;
data.periodLabels    = periodLabels;
data.geoLabels       = geoLabels0;
data.surfaces        = surfaces;
data.population      = population;
data.incomingCommuters = incomingCommuters;
data.touristBeds       = touristBeds;

% Save to output file
save(runDataFile, 'data', '-v7.3');

% Save things that are useful for dynamic prediction and counterfactuals
base.Y             = Y; % (sparse) T x K1 x K2
base.Xdynamic_staticInputs = Xdynamic_staticInputs; % cell(NumParts,1)
base.laggedConsumersOrig = reshape(laggedConsumersOrig, [T K1]);  % T x K1
base.laggedConsumersDest = reshape(laggedConsumersDest, [T K2]);  % T x K2
base.laggedProviders  = reshape(laggedProviders,  [T K2]); % T x K2
base.popAtRisk     = popAtRisk; % T x K1
save(dynPredInputFile, 'base', '-v7.3');

clear;
end